# 1
library(Sleuth2)
attach(case0902)
names(case0902)
library(MASS)

lbrain = log(Brain)
lbody = log(Body)
lgest = log(Gestation)
llitter = log(Litter); detach()

#X = cbind(rep(1,96),lbody,lgest,llitter)
#Vbeta = solve(t(X)%*%X)
#bhat = Vbeta%*%t(X)%*%lbrain
#
#s2 = as.numeric(t(lbrain-X%*%bhat)%*%(lbrain-X%*%bhat))/92
#sigma2 = 92*s2/rchisq(1000,92)
#betas = matrix(0,1000,4)
#for (i in 1:1000) {
#	betas[i,] = as.vector(mvrnorm(1,as.vector(bhat),
#	(Vbeta*sigma2[i])))
#}
#
#apply(betas,2,mean)
#sqrt(apply(betas,2,var))
#apply(betas,2,quantile,c(0.025,0.975))

m1 = lm(lbrain ~ lbody)
m2 = lm(lbrain ~ lgest)
m3 = lm(lbrain ~ llitter)
m4 = lm(lbrain ~ lbody + lgest)
m5 = lm(lbrain ~ lbody + llitter)
m6 = lm(lbrain ~ lgest + llitter)
m7 = lm(lbrain ~ lbody + lgest + llitter)

bic1 = BIC(m1)
bic2 = BIC(m2)
bic3 = BIC(m3)
bic4 = BIC(m4)
bic5 = BIC(m5)
bic6 = BIC(m6)
bic7 = BIC(m7)

bicv = cbind(bic1,bic2,bic3,bic4,bic5,bic6,bic7)
postModel = 1/7 * exp(-bicv) / sum(1/7 * exp(-bicv))

# 2
countBicy = c(16, 9, 10, 13, 19, 12, 1, 4, 9)
countOther = c(58, 90, 48, 57, 103, 12, 1, 4, 9)
N = countBicy + countOther

postTheta = matrix(0,9,100000)
for (i in 1:9)
{
    postTheta[i,] = rbeta(rep(1,100000), 1+countBicy[i], 1+N[i]-countBicy[i])
}
phiY = 1/5 * colSums(postTheta[1:5,])
phiN = 1/4 * colSums(postTheta[6:9,])
meanDiff = mean(phiY - phiN)
PostQuan = quantile(phiY - phiN, c(0.025, 0.975))
